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A statistical model of fragmentation of aggregates is proposed, based on the stochastic propa- 
gation of cracks through the body. The propagation rules are formulated on a lattice and mimic 
two important features of the process - a crack moves against the stress gradient and its energy 
depletes as it grows. We perform numerical simulations of the model for two-dimensional lattice 
and reveal that the mass distribution for small and intermediate-size fragments obeys a power-law, 
F(m) oc m -3 / 2 , in agreement with experimental observations. We develop an analytical theory 
which explains the detected power-law and demonstrate that the overall fragment mass distribution 
in our model agrees qualitatively with that, observed in experiments. 
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Introduction. Fragmentation processes are ubiquitous 
in nature and play an important role in many industrial 
processes. Numerous examples range from manufactur- 
ing comminution to collision of cosmic bodies in space. 
Hence, much effort has been devoted to comprehend the 
nature of fragmentation in different branches of research 
- geophysics Jl], Hi ; astrophysics SB, engineering [j| and 
material- [(| or military science [3,111 ■ 

An intriguing common property of fragmentation is 
a power law mass-distribution of the fragments which 
seems to be independent of the spatial scale or nature of 
the parent bodies: When heavy ions collide with a target, 
or asteroids suffer a high-speed impact, both produce a 
power law size (or charge) distribution of debris. This 
law has been reported in numerous experimental studies, 
e.g. [Mil. 

The theoretical description of fragmentation is devel- 
oping along two different lines: one, based on continuum 
mechanics, another - on the general statistical methods, 
e.g. HE El. Finite element analysis (FEM) [13|, nu- 
merical simulations of agglomerate collisions with smooth 
particle hydrodynamics (SMH) [2l|, or discrete element 
method (DEM) 0, Il7| — these are the current tools to 
treat continuum mechanical problem of colliding mate- 
rial bodies. Alternatively, statistical methods, as random 
wal k |22l| or stochastic simulation of crack-tip trajecto- 
ries [23| were used to model fractal crack formation and 
propagation. Furthermore, there were attempts to tackle 
the problem analytically. In particular, a semi-empirical 
approach to the physics of catastrophic breakup of cos- 
mic bodies has been developed 
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similarly, a micro- 
structural approach has been adopted [25[ to model the 
fracture generation and propagation in concrete. 

The observed universality of the fragmentation law for 
the objects, drastically different in material properties 
and dimension, most probably implies a common physi- 
cal principle inherent to all fragmentation processes |26[ . 
Moreover, it is reasonable to assume that this principle 



is of statistical nature. This motivates the analysis of a 
model, which being very simple on the microscopic level, 
reflects the most prominent features of the process and 
adequately represents its statistical properties. In the 
present Letter we propose a novel statistical model for 
fragmentation of aggregates - macroscopic bodies, com- 
prised of a large amount of smaller macroscopic con- 
stituents. The problem of aggregate fragmentation arises 
in many areas of science and technology, in particular in 
planetary science, where the availability of an adequate 
fragmentation model is crucial for understandin g of the 
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formation and evolution of planetary rings, e.g. 
We show that our model, based on a few very simple 
physical rules reproduces the main aspects of the frag- 
mentation process and gives the power-law distribution 
for the debris size. Moreover, it qualitatively agrees with 
the experimental data. We perform numerical simula- 
tions and develop a simple analytical theory, which ex- 
plains the observed power-law distribution for the frag- 
ments. 

Formulation of the model. We consider aggregates, 
which are not too large to neglect the gravitation forces 
and assume that particles are kept together by the ad- 
hesive bonds. These are significantly weaker than the 
forces attributed to chemical bonds, hence the formation 
of cracks occurs only along the contact lines between the 
particles. We assume that all particles are of the same 
size and material. Then to break any adhesive contact, 
the amount of energy, equal to E\>, is required; to de- 
stroy simultaneously n bonds one needs ??-timcs larger 
energy, nE^. The quantity E\, depends on the particles 
size, their surface tension and material parameters, e.g. 
[30l | . To simplify further the problem, we consider its 
two-dimensional version and assume that the spherical 
constituents form a square lattice, Fig. [T] 

In a collision between aggregates, which causes their 
fragmentation, the energy of the relative motion -E C oii 
is transformed into the surface energy of broken bonds 
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and the kinetic energy of debris. We will ignore the lat- 
ter part of the energy and explicitly consider the for- 
mer one. Therefore, the problem of fragmentation in our 
model is, essentially, reduced to the problem of distri- 
bution of the energy E co i\ among the broken adhesive 
bonds. We assume the bonds being destroyed due the 
crack's propagation. Namely, we adopt following rules 
for the crack growth: (i) A crack originates on the sur- 
face, on the impact site, where the maximum stress is 
expected. It propagates (in average) away from the sur- 
face against the stress gradient, that is, "downhill" the 
load, and never "uphill", (ii) Each time-step the crack 
elongates by a single neighboring bond, consuming the 
energy Ey,, while the crack tip performs a random walk 
with the direction randomly chosen among two or three 
possible ones, as explained in Fig. [TJ (iii) The propaga- 
tion of a certain crack terminates if its tip arrives at the 
surface, or meets another crack, mimicking a crack bi- 
furcation. When a crack terminates and the rest-energy 
allows for further bond-breaking, a new crack is initi- 
ated randomly on the surface of the aggregate, or bifur- 
cates from another crack, (iv) A mechanism limiting a 
total crack energy is introduced. This is because crack 
generation and propagation needs a sufficient load gradi- 
ent which decreases with distance from the impact site. 
Thus, we assume that if the energy assigned to a crack, 
chosen randomly between zero and a maximum E cross , 
is exhausted, the crack stops and another one is created 
instead. We called this "crossing factor", (v) The break- 
ing process continues as long as the remaining impact en- 
ergy, -finest = E co n — nEh (n - current number of broken 
bonds) is sufficient to break further bonds. It terminates 
whenever the total impact energy E co \\ is dissipated in 
cracks. 
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FIG. 1: (Color on line). A model aggregate is composed of 
I x I identical particles on a square lattice. To break a bond 
between particles the energy E^ is needed. A crack propagates 
either away from the loaded surface (vertically-down on the 
plots), against the stress gradient, or laterally, but never along 
the gradient. All possible directions have equal probability. 
In the case (a) there exist three allowed directions, each with 
the probability ps = 1/3: one away from the surface and two 
lateral directions. In the cases (b) and (c) only two directions 
with the probability p2 = 1/2 are possible, since crack can 
not move back or "uphill" the load. 

In Fig. [2] a typical fragmentation pattern, obtained 
with the use of the above fragmentation rules, is shown. 
Next, we analyze the statistical properties of debris. 

Fragment mass distribution. We define a fragment as 
a collection of constituent particles connected by adhc- 





FIG. 2: Left panel: A typical fragmentation pattern. The ag- 
gregate, composed of one million particles on the 1 000 x 1 000 
lattice, is broken in 3958 pieces. The total breakage energy 
is .Ecoii = 82 010£'b, the crossing-factor is E CIOBS = 3 000 Eb- 
Right panel: The fragmentation model may be mapped on 
the one-dimensional model of diffusing and annihilating par- 
ticles. The areas under the particle's trajectories correspond 
to the fragment areas of the primary fragmentation model. 



sive bonds where a single particle consitutes a mass unit. 
Hence, for a fragment of m particles we assign the mass 
to. We have performed numerical experiments, control- 
ling the size of the aggregate, the total energy spent in 
the process and the crossing factor. To improve the frag- 
ments' statistics we perform many runs with the same 
set of parameters. 

Fig.[3]shows the fragment mass-distribution. The most 
part of the distribution, except for very large fragments, 
may be accurately fitted by the power-law 



F(m) oc to a 



(1) 



Interestingly, the exponent a ~ 1.5 is the same as in 2D 



egg-shell crushing experiments [14 
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FIG. 3: (Color online). Fragment mass distribution for a 
100 x 100 lattice matrix with total energy of 4 000 Eb and a 
crossing factor of 200 -Eb- The data have been collected for 
10 000 runs. The line represents a power law mass-distribution 
with an exponent ~ —3/2. The inset shows the distribution 
function P(x, t) oc i _3//2 , which gives the probability that two 
particles, initially separated by distance x, meet at time t; it 
illustrates the affinity of the two models (see the text). Points 
- numerical data for x = 20, line - theory. 



Analytical model. The observed power-law fragment 
size distribution may be obtained analytically, if we 
notice that the proposed fragmentation model may be 
mapped onto the one-dimensional model of diffusing and 
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annihilating particles, where the particles correspond to 
the tips of the cracks. Indeed, the direction normal to 
the surface, that is, the direction of the most rapid de- 
cay of the stress, say axis y (vertical downward in Fig. 
[T]), may be mapped onto the time axis, since the reverse 
motion along this axis is forbidden. Then the location of 
the crack tips along the lateral direction, say along axis 
x (horizontal in Fig. [TJ corresponds to the location of 
fictitious "particles" on the line. One step along y axis 
corresponds to one time step, while one step along x axis 
corresponds to a "particle" jump on the line (without loss 
of generality we can assume unit time and space steps). 

It is easy to see that the fragmentation model with the 
rules illustrated in Fig. Q] is equivalent to the following 
model: (i) Each particle on a line at each time step can 
remain at the same site with the probability pq = 1/3, 
move one site left or right with the probability p\ = 1/6, 
two sites left or right with p2 = 1/12, etc., k sites left or 
right with the probability pk = 1/(3 • 2 k ); note that po + 
2j2kPk = 1- (ii) When two particles meet, one of them 
annihilates while the other one continues to move with 
the same rules, (iii) The fragment sizes of the primary 
fragmentation model is equal to the area confined by the 
trajectories of the particles and the axis x, Fig. [5J 

First we show that the "particles" perform one- 
dimensional diffusion motion. Indeed, due to the lack 
of memory, each time step does not depend on the pre- 
vious one (note that the original problem does have a 
memory with respect to previous steps). Therefore, the 
mean square displacement < [Aa;(<)] > is a sum of these 
quantities for each step, (A 2 ^), where 



~ 2* 2 
^ 3 • 2 l ~ 



1=0 



1=0 



that is, {[Ax(t)] j = At = 2Dt. Hence, we have a diffu- 
sive motion with the diffusion coefficient D = 2. 

Now we compute the distribution of the fragments' ar- 
eas. These correspond to the areas of the figures on the 
x — t plane confined by the x-axis and the diffusive trajec- 
tories of pairs of particles, which start to move at t = 0, 
and meet at t > 0; the initial separation of the particles 
is x. The area of the figures, that is, the fragment mass 
scales as m cx tx. Hence, the mass distribution reads: 



F(m) 



dx / dtS(m — jtx)xe >cx P(x,t). 



(2) 



where the coefficient 7 accounts for the surface density 
and a geometric factor (1/2 for triangles), which relates 
a fragment area and its dimensions x and t. The coef- 
ficient h is the line density of particles on the x axis at 
the starting time t = (i.e. the density of crack tips 
which start at y = 0); the random initial distribution of 
particles (crack tips) implies the Poisson distribution for 
the initial inter-particle distance, x exp(— xx). P(x,t) 



gives the probability of two diffusing particles, initially 
separated by distance x, to meet at time t\ where one 
of the particles annihilates. P(x, t) may be written as 
P{x,t) = jj-PsurvOM), where 



f°° 1 / 

^ t)= Jovmi^ 



-{y-xY/SDt 



-(y+x) 2 /8Dt 



dy 
(3) 

is the survival probability, i.e. that both particles have 
not annihilated before time t, provided they started to 
diffuse at t = separated by the distance x. It is ob- 
tained from the solution of the diffusion problem with 
the adsorbing boundary condition for the relative mo- 
tion of two particles, with the double diffusion coefficient 
2D (see e.g. [1H). Indeed the integrand in Eq. © gives 
the probability distribution of the inter-particle distance 
y at time t, if the initial separation was x, provided this 
probability is identically zero for y = (the annihila- 
tion condition). Integrating this probability over all dis- 
tances y > gives the survival probability; differentiation 
of P SUTV (x,t) with respect to time yields the probabil- 
ity that the particles meet exactly at time t. Substitut- 
ing P(x,t) = x exp(-x 2 /8Dt) /V32Di 3 , obtained from 
Eq. ([3]) into Eq. ([2]), we arrive at the expansion, 



F{m) = A m- 3/2 + Aim _5/a + 



-5/2 



(4) 



where A Q = 7 1 / 2 r(5/2)/8xr 3/2 and generally, A t = 
(-l)'r(31 + 5/2)(7/^) 3 ' +1 /2/(8x)(4 7 ) 2/ Z! is the coeffi- 
cient at to - ( 3 / 2+ '). Hence, Eq. (fj| explains the power-law 
distribution (fTJ) obtained numerically. The above expan- 
sion holds, however, not for very small masses, to keep 
the continuum diffusion approach valid, and not for very 
large masses, to ignore the effects of energy depletion and 
finite size effects of the fragmenting pattern. 

Comparison with experiments. So far we considered 
the size distribution of fragments which mass is much 
smaller than the initial mass of the parent body. In 
this limit the numerically detected and explained power- 
law distribution coincides with the experimental one 14 1 . 
Surprisingly, our model can qualitatively and sometimes 
even semi-quantitativcly reproduce the whole size distri- 
bution of fragments in the impact experiments. This is 
illustrated in Fig.|4]for the experiments of exploding egg- 
shells [lj] and in Fig. [5] for the colliding basalt spheres 
[l(]| . For the former case our model demonstrates an over- 
all qualitatively correct behavior for F(m), while for the 
latter case one can achieve a quantitative agreement with 
the most part of the experimental fragment size distribu- 
tion. It is worth noting that the agreement is obtained 
without any fitting or scaling of our data. Moreover, the 
fact that our simple two-dimensional model can describe 
the fragmentation statistics for three-dimensional bodies 
indicates, that the main ingredients of our model - the 
diffusive propagation of cracks against the stress gradient 
and the depletion of energy with the cracks' growth - are 
indeed the basic features of a fragmentation process. 
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FIG. 4: (Color online). Results of our model (red circles) 
compared with the experimental data (green squares) of Her- 
rmann et al. [3 • In simulations we used a 600 x 600 matrix, 
the total energy of -E co ii = 40 000 Eb, the crossing-factor of 
E crosa — 3 000 Eb and the number of runs is 200. Note the 
lack of any fitting or scaling of our data. 
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FIG. 5: (Color online). Results of our model (red circles) com- 
pared with the experimental data (green circles) of Nakamura 
and Fujiwara In simulations we used a 80 X 80 matrix, 
the total energy of E co \\ — 1500 Eb, the crossing-factor of 
E CTOaa — 100 Eb\ the number of runs is 10. 



In conclusion we suggest a simple fragmentation 
model, based on a lattice random walk of a crack with 
propagation rules which mimic a real fragmentation pro- 
cess: A crack moves against (or laterally to) the stress 
gradient and the energy of the crack exhausts along 
its growth. We performed numerical simulations of the 
model and observed that the mass distribution of small 
and intermediate fragments obeys a power-law. The ex- 
ponent of the power law, equal to —3/2, agrees with the 
reported experimental value. We develop an analytical 
theory which explains the nature of the power-law in the 
fragmentation statistics and demonstrate that the pro- 
posed model gives a qualitative description for the overall 
fragment size distribution observed in experiments. 
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